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Abstract: 

We study a recent model for calcium signal transduction. This model displays 
spiking, bursting and chaotic oscillations in accordance with experimental re- 
sults. We calculate bifurcation diagrams and study the bursting behaviour 
in detail. This behaviour is classified according to the dynamics of separated 
slow and fast subsystems. It is shown to be of the Fold-Hopf type, a type 
which was previously only described in the context of neuronal systems, but 
not in the context of signal transduction in the cell. 



1 Introduction 

Nonlinear dynamics has attracted much attention in several scientific fields, 
predominantly in physics. However, oscillations and chaos have also been 
discovered within several levels of abstraction in biology and (bio) chemistry. 
These levels range from populations of organisms (predator and prey systems) 
[T], individual organisms (circadian and ultradian clocks) [2] down to bio- 
chemical reaction pathways (glycolysis) 3j and even individual (bio) chemical 
reactions (e.g. PO reaction) [I]. In most of these cases, chaotic oscillations 
have been observed in addition to periodic behaviour. However, the route to 
chaos can be only studied in experimentally well tractable systems. 

In addition to simple periodic and chaotic oscillations, complex periodic 
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and quasiperiodic oscillations have been observed. Among the complex os- 
cillations, bursting behaviour has been studied thoroughly in the biological 
context, since it was observed in important processes such as nerve signal 
conduction [5] and signal transduction [H] within the cell. Several types of 
bursting have been characterized. 

Here, we study a recent model for calcium signal transduction [7 j. Cal- 
cium ions play a central role as second messengers within the cell. Following 
the stimulation of agonist receptors at the cell membrane, a cascade of events 
is triggered which finally leads to the massive liberation of calcium from in- 
tracellular stores. Calcium in turn influences multiple processes like gene 
expression, diverse enzymatic reactions, vesicular transport and embryogen- 
esis [Hill. 

In agreement with experimental results, the model displays simple peri- 
odic (spiking), bursting and chaotic oscillations. We compare the bursting 
behaviour of the complete model to an independent dynamics of its slow and 
fast subsystems. We study the bifurcation behaviour of the model with one 
slow variable being used as the bifurcation parameter, considering all other 
variables to be slaved. This allows the classification of bursting, an approach 
first used by J. Rinzel ^U] for bursting electrical activity in models of nerve 
cells. We classify the bursting behaviour found in the model under investi- 
gation to be of the Fold-Hopf type according to the comprehensive scheme 
proposed in [TT| . So far, this type of bursting had been described in the 
context of neural activity. Here, we report the first observation in a model 
for calcium signal transduction. Reasons for this behaviour in terms of the 
connectivity of the system are discussed. 

2 Models for Calcium Dynamics 

Many calcium signal transduction models have been developed (for a review 
sec [12J). The one studied here and illustrated in Fig. ^describes the following 
scenario: Upon binding of the agonist, the G a -subunit of the receptor coupled 
G-protein is activated. Variable a denotes the concentration of activated G a - 
subunit. This in turn activates phospholipase C. The concentration of the 
activated form PLC* is denoted by h. The activated G a -subunit and PLC* 
cause the liberation of calcium from the endoplasmic reticulum (ER) and 
from the extracellular space. The calcium concentrations in the cytosol and 
in the endoplasmic reticulum are described by the variables c and d. Calcium 
is pumped back to the ER and to the extracellular space. Moreover, different 
mechanisms lead to the inactivation of a and b over time. 

Variables b and c exhibit negative feedback on variable a by acceleration of 
the inactivation processes. In addition, variables a and c show autocatalytic 
behaviour [7j. The autocatalysis in c is an integral compound of a term that 
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Figure 1: Schematic drawing of the cell with a membrane bound receptor 
(top) and the endoplasmic reticulum (ER). Dotted arrows denote the mod- 
elled influences on the Ca 2+ fluxes (solid arrows). 



describes the activity of a calcium channel at the endoplasmic reticulum 
(ER). We modified this term compared to the original model in order to 
make it more realistic. The change mainly ensures the reversibility of calcium 
diffusion through the channel. The qualitative dynamics of the model is not 
affected by this change. 

The model thus has the following form: 



da k 3 ab k$ac . , 

— = h + k 2 a — — (1) 

dt a + ki a + kq 

^ = kra-^T (2) 

dt b + k 9 w 

dc k w cb 4 (d - c) k 1A c k 16 c 

dt = WTWi +kl2b + kl3a -^Th 5 -^Th 7 (3) 

dd k w cb A (d — c) k^c . . 

dt = 6 4 + k\ x + c + k 17 ^ 

This model exhibits spiking and bursting, as well as chaotic oscillations 
in dependence on agonist stimulation that is proportional to the parameter 
k 2 [7j. Moreover, it can be reduced to a 3 variable model which still captures 
the essential dynamics of the 4 variable model. This reduced model has the 
following form: 



da , k^ab k$ac 

— = k x + k 2 a — — 

dt a + K4 a + fc 6 



(5) 
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Again, variables b and c exhibit negative feedback on a. The autocatalysis 
in a is still present. However, the autocatalysis in c is not necessary. 



In order to perform a bifurcation and linear stability analysis, the numerical 
continuation software AUT097 ^3] is employed. AUT097 computes fixed 
points and limit cycles of sets of ordinary differential equations. It calculates 
branches of these solutions under parameter variation, detects bifurcation 
points and allows to switch to a new branch emerging at a bifurcation. It 
also monitors eigenvalues respectively Floquet multipliers of the solutions to 
obtain their linear stability. This method is used for both the analysis of the 
complete system as well as for its separated slow and fast subsystems. 

The classification of the bursting behaviour is carried out according to 
the following procedure: The bursting behaviour consists of two consecutive 
phases, a silent phase with slow changes of the variables and an active phase 
with rapid oscillations of some of the variables whereas others remain slowly 
varying. Hence, during the active phase, the individual variables change on 
different time scales and may be grouped into a slow and a fast subsystem, 
respectively [Ej. The encountered bifurcations of the fast subsystem that 
initiate and terminate the active phase classify the type of bursting according 
to PI]. 

4 Results 

First, we will focus on the 3 variable model ©-(0) and later compare its 
bifurcation structure to the 4 variable model. We choose k 2 as the bifurcation 
parameter. It denotes the strength of agonist stimulation of the receptor and 
it is experimentally well studied. All other parameters are kept fixed. The 
bifurcation diagrams in Fig.|2]show the maximum and minimum of c(t) versus 
the parameter k 2 . 

At k 2 < 1-3, the steady state is stable and c increases monotonically with 
the stimulation k 2 - At k 2 = 1.3, a supercritical Hopf bifurcation renders 
the steady state unstable and stable periodic oscillations coexist with the 
unstable steady state for k 2 > 1.3. As k 2 is increased in the interval 1.3 < 
k 2 < 2.92 the oscillations become increasingly complex and develop secondary 
maxima in c(t). In Fig. ^examples of spiking (a)-(c) at k 2 = 2.0 and bursting 
(d)-(f) behaviour at k 2 = 2.9 are compared. 



3 Methods 
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Figure 2: Bifurcation diagrams of complex calcium oscillations (thick curves 
show maximum and minimum c(t)) and steady state behaviour (thin curves 
in (a,b)) in Eqs. ©-© as function of the agonist level k 2 . Solid (dashed) 
curves give stable (unstable) solutions. Symbols denote Hopf (open square), 
saddle-node of periodic orbits (filled triangle) and period doubling (open 
diamond) bifurcations, (a) covers the whole parameter range of interest 
whereas (b) is a close-up of the steady state (thin) and the maximum c 
(thick curve) in the parameter range where periodic bursting is unstable and 
chaotic behaviour has been observed, e.g. at k 2 = 2.9259 in Fig.7 of 7J. (c) 
shows one of the period doubling cascades with branches of 1 (thick), 2, 4 and 
8 irregular bursts (thin) per period. See the text for details. Parameters k\ = 
0.212, k 3 = 1.52, fc 4 = 0.19, k 5 = 4.88, k % = 1.18, k 7 = 1.24 } k 8 = 32.24, k 9 = 
29.09, ki 3 = 13.58, ku = 153, ki 5 = 0.16 are the same as in Fig.7 of jjj. 
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At k 2 = 2.9258, the branch of stable bursting undergoes a period doubling 
bifurcation and turns unstable. Stable oscillations are recovered after a se- 
quence of saddle-node and period doubling bifurcations of periodic orbits as 
shown in Fig. |2Jb). Each hysteresis reduces the number of secondary spikes 
per burst by one. The intervals of increasing maximum of c(t) with k 2 are 
unstable due to the saddle-node bifurcations. Those of decreasing maximum 
are destabilized by pairs of period doubling cascades that occur very close to 
the saddle-node bifurcations (compare panels (b) and (c) in Fig. |2J). 

Fig. |2fc) gives a close-up of a period doubling cascade from Fig. I2fb) that 
is the reason for chaotic behaviour within the parameter intervals between 
pairs of such cascades. We only show 4 period doubling bifurcations with 
branches of 1 (thick), 2, 4 and 8 bursts per period. Due to increasing com- 
putational effort, the emerging branch with 16 irregular bursts per period is 
not computed further. This chaotic behaviour has been reported earlier at 
k 2 = 2.9259 (Fig. 7 in [7j). 

After the periodic oscillations changed back to simple spiking they termi- 
nate in the second supercritical Hopf bifurcation at k 2 = 3.0. For k 2 > 3.0, 
the steady state is stable and corresponds to over-stimulation with constant 
high calcium concentration in the cytosol. At large k 2 , the model displays a 
third Hopf bifurcation where again stable periodic oscillations emerge. How- 
ever, this feature lies outside the parameter range of interest. 

From Fig.E^d)-(f) it is evident that the smooth oscillations of b(t) evolve 
on a slower time scale than the secondary spikes of a(t) and c(t). The same 
solution is given by the dotted curve (with circles) in Fig.0J Here, the dynam- 
ics in the three-dimensional phase space (a, b, c) is displayed by projection 
onto the (b, a) and (6, c) planes. 

In the following, the slow subsystem Eq. © and the fast subsystem 
Eqs. ©,(0) are studied separately. Setting Eq. © equal to zero one de- 
rives the nullcline 

k 8 /k 7 b 

ao{b) = bTh (8) 

where db/dt > for a > ao(b) and db/dt < for a < ao(b) as follows from the 
first term on the right hand side of Eq. © • 

The bifurcation diagrams of the fast subsystem Eqs. ©,© with b as 
bifurcation parameter are superimposed in Fig. 0] As b is varied, the steady 
state solution (thin curves) for a and c shows a hysteresis limited by two 
saddle-node bifurcations (filled triangles). The lower branch is stable and 
the middle branch unstable. On the upper branch, a Hopf bifurcation (open 
square) gives raise to fast oscillations. In the present model, the time scale 
separation is finite and the solution of the full system ©-© slightly deviates 
from the curves derived from the fast subsystem. 

We gain a qualitative understanding of the complex structure of burst- 
ing behaviour if we follow the circles in clockwise direction (arrow) as time 
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Figure 3: Examples of spiking (a)-(c) at k 2 = 2.0 and bursting (d)-(f) at 
k 2 = 2.9 in Eqs (JHJ-fjTJ). The ordinate labels apply to both panels in each 
row. The bursts (d)-(f) consist of a silent phase with low values of a(t) and 
c(t) while b(t) decreases and an active phase with rapid oscillations of a(t) 
and c(t) when b(t) monotonously increases. The dynamics of a(t) and c(t) 
occur on a fast and the dynamics of b(t) on a slow time scale. Parameters 
are the same as in Fig. |21 
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progresses. Starting in the silent phase of the burst where a and c are low 
(compare Fig. E^d)-(f)) the dynamics are close to the lower stable branch of 
a(b) and c(b). db/dt < holds, since the nullcline (thick dashed curve) lies 
above. Hence, b will slowly decrease and a and c will adiabatically follow the 
lower branch in Fig. 0] until the lower saddle-node bifurcation is exceeded. 
Then, no stable steady state exists any more and the dynamics approach the 
fast oscillation (thick curve) at larger values of a and c. Now a > a and b 
increases. 

The crossing of the saddle-node bifurcation triggers the active phase of 
the burst which starts with a large spike in a(t),c(t) followed by secondary 
spikes of decreasing amplitude. The spikes correspond to the branch of fast 
oscillations in Fig. 0] At larger b these oscillations become smaller and con- 
tinuously vanish in the supercritical Hopf bifurcation. However, since the 
time scale separation is finite, the complete dynamics passes the Hopf bifur- 
cation and returns via damped oscillations to the stable focus on the upper 
branch. This feature has been analysed in detail in a model of enzyme kinet- 
ics that shows the same type of bursting [14J. After a short plateau (upper 
stable steady state) the saddle-node bifurcation at large b is exceeded. The 
lower steady state is approached and b starts to decrease. This completes 
one period of the burst. 

Since the active phase starts at a saddle-node bifurcation (fold) and ends 
on a branch of stable foci provided by a supercritical Hopf bifurcation, the 
present dynamics is an example of Fold-Hopf bursting following the classifi- 
cation by E. M. Izhikevich [TT] . 

We have studied the 4 variable model ©-(JU) in the same way as the 
3 variable model (J5])-({7j). The two models qualitatively show the same be- 
haviour. In Fig. we show the bifurcation diagram of the 4 variable model. 
The solutions have a similar shape and the bifurcation diagrams of the fast 
subsystems are qualitatively the same. In the two models b(t) is the only 
slow variable and the remaining variables constitute the fast subsystem. The 
type of bursting is the same in the two models. 

5 Discussion 

We studied the bifurcation and bursting behaviour of a recently proposed 
model for calcium signal transduction and classified the bursting behaviour to 
be of Fold-Hopf type. This Fold-Hopf type of bursting has also been observed 
in the electrical activity of pyramidal cells of the cat hippocampus ^5] and 
in models of the bursting electrical activity in pancreatic /3-cells ^H]- The 
authors of the latter paper call it "tapered" bursting. The models have been 
reduced to the Lienard form which has been studied intensively [T71 ll8|. The 
same type of bursting was reported and called "Type V" [T7] . However, the 
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Figure 4: Solutions of the fast subsystem (a(t,b),c(t,b)) at k 2 = 2.9 as in 
Fig. 0(d)- (f), projected onto the (b, a) and the (6, c) plane. Thin (thick) 
curves stand for steady states (maxima and minima of fast oscillations) and 
dashed (full) curves indicate unstable (stable) solutions. The dotted curve 
(with circles) represents the projection of the solution of the full system 
©-(0) which is parametrized by time in the clockwise direction following 
the arrow. Small irregularities at low a and b reflect the finite numerical 
accuracy. The thick dashed curve in (a) represents the nullcline b t = 
of the slow dynamics. b(t) increases (decreases) above (below) this curve. 
Parameters are the same as in Fig. [21 
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Figure 5: The bifurcation diagram of the 4 variable model CD)-® is very 
similar to that of the 3 variable model ©-©• The branch of unstable 
oscillations (thick dashed) has the same structure as in Fig. El . Curves have 
the same interpretation as in Fig. |2l and parameters are k\ = 0.09, = 
0.64, k 4 = 0.19, k 5 = 4.88, k 6 = 1.18, k 7 = 2.08, k 8 = 32.24, k 9 = 29.09, k 10 = 
5.0, k u = 2.67, k 12 = 0.7, k 13 = 13.58, k 14 = 153, k 15 = 0.16, k 16 = 4.85, k 17 = 
0.05 . 



classification scheme for bursting by E. M. Izhikevich [11 follows a bottom- up 
approach and provides a systematic nomenclature. It is therefore preferred 
by us. 

In the purely biochemical context, the only previous example for Fold- 
Hopf type bursting is an abstract enzymatic system [T^. This system is 
also subject of the study mentioned above ^3]. It can be described by three 
variables and has the following form: 



a' = K-5(f)(a,(3) (9) 
& = q l 8<j>{a^)-8v{f3, 1 ) (10) 
i = q 2 5u((3, 1 )-k sl (11) 

It is interesting to compare the topology of this abstract biochemical 
system with the topology of the calcium signal transduction system studied 
here. The abstract system involves two negative feedback loops of variable 
(3 on a as well as 7 on (3. It also contains two autocatalysis of the variables 
7 and (3. 

The 4 variable model studied here also contains two negative feedback 
loops as well as two autocatalytic species. However, the 3 variable model 
which displays the same dynamics contains only one autocatalytic species. 
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Therefore, the overall topology of the two systems is similar but not the 
same. 

In [211], the model (|§J)-(fTT|) has been extended by two quickly relaxing 
variables that account for the finite time scale of allosteric transitions in the 
enzymatic system. The bursting behaviour was shown to depend on the 
chosen separation of time scales with the original dynamics flU|)- (jll|) being 
conserved in the limit of fast allosteric transitions. Our 3 variable and 4 
variable models also assume fast allosteric transitions and so far other exper- 
imental indications are not available. On the other hand, our models possess 
a different topology and the results for slow allosteric transitions [20] do not 
immediately apply to the present situation. 

In general, negative feedback loops and simple autocatalysis are very 
common in biochemical systems. Cooperative or allosteric behaviour in en- 
zymatic systems can be considered to be an autocatalysis for certain ranges 
of parameters. This behaviour is involved in metabolic regulation and sig- 
nal transduction. It offers the possibility to fine tune biochemical systems, 
but also to display a wealth of different dynamic properties. Since the basic 
requirements for this behaviour are so abundant in the biochemical network 
in living cells, it is likely that there are many more examples which have not 
been discovered yet. Bursting behaviour in biology seems to be an integral 
part of the information processing machinery. 

How information is processed in the cell is still not very well understood. 
However, it has been studied in some detail in the context of calcium signal 
transduction. It has been shown that the frequency of calcium oscillations 
encodes information [2B 1221 12H1 121] • Fundamental differences in the topology 
of the signal transduction pathway also offer the possibility to give rise to 
fundamentally different dynamical behaviour and encode information in this 
way jjj. In this context, it is of importance which type of bursting behaviour 
is displayed. Upon variation of the model parameters the locations of bifur- 
cation points shift along the steady state branches in the fast subsystem. For 
different types of bursting these quantitative changes may trigger different 
qualitative changes of the bursting behaviour. 

Other models of calcium dynamics have recently been studied that also 
show bursting behaviour but of Fold-Sub- Hopf hysteresis [22] and Sub-Hopf- 
Fold-Cycle type [2E], respectively. In these models the amplitude of the 
secondary spikes increases as the end of the active phase is approached. This 
corresponds to unstable foci in their fast subsystems that are provided by 
subcritical Hopf bifurcations. With our model, the secondary spikes decrease 
in amplitude. More detailed experimental studies in the future will reveal 
which type(s) of bursting are predominant in calcium signal transduction. 

Acknowledgements: 

UK and MO thank the Klaus Tschira Foundation for funding. 



11 



References 

[1] M. J. Begon, L. Harper and C. R. Townsend, Ecology: Individuals, 
Populations, and Communities, 3rd edition, Blackwell Science Ltd. 
Cambridge, MA (1996). 

[2] L. N. Edmunds (ed.) Cell Cycle Clocks, Marcel Dekker, New York 
(1984). 

[3] A. Ghosh and B. Chance, Biochem. Biophys. Res. Commun. 16 (1964) 
174. 

[4] L. F. Olsen and H. Degn, Nature 267 (1977) 177. 

[5] X. J. Wang and J. Rinzel, Brain Theory and Neural Networks (M. A. 
Arbib, ed.), The MIT press, Cambridge, MA (1995). 

[6] C. J. Dixon, N. M. Woods, T. E. Webb and A. K. Green, Biochem. J. 
269 (1990) 499. 

[7] U. Kummer, L. F. Olsen, C. J. Dixon, A. K. Green, E. Bornberg-Bauer 
and G. Baier, Biophys. J. 79 (2000) 1188. 

[8] M. J. Berridge, Nature 361 (1993) 315. 

[9] M. J. Berridge, M. D. Bootman and P. Lipp, Nature 395 (1998) 645. 

[10] J. Rinzel, in Mathematical topics in population biology, morphogene- 
sis, and neurosciences, Eds. E. Teramoto and M. Yamaguti, Springer, 
Berlin, (1987). 

[11] E. M. Izhikevich, Int. J. Bifurcat. Chaos 10 (2000) 1171. 

[12] J. Sneyd, J. Keizer and M. J. Sanderson, FASEB J. 9 (1995) 1463. 

[13] E. Doedel, A. Champneys, T. Fairgrieve, Y. Kusntsov, et al, 
AUT '097 '-.Continuation and bifurcation software for ordinary differ- 
ential equations (Concordia University, Montreal, 1997). 

[14] L. Holden and T. Erneux, SIAM J. Appl. Math. 53 (1993) 1045. 

[15] E. R. Kandel and W. A. Spencer, J. Neurophysiol. 24 (1961) 243. 

[16] P. Smolen, D. Terman and J. Rinzel, SIAM J. Appl. Math. 53 (1993) 
861. 

[17] G. de Vries, J. Nonlinear Sci. 8 (1998) 281. 

[18] M. Pernarowski, SIAM J. Appl. Math. 54 (1994) 814. 



12 



[19] O. Decroly and A. Goldbeter, J. Theoret. Biol. 124 (1987) 219. 

[20] M. Kaern and A. Hunding, J. Theoret. Biol. 198 (1999) 269. 

[21] W. Li, J. Llopis, M. Whitney, G. Zlokarnik and R. Y. Tsien, Nature 
392 (1998) 936. 

[22] R. E. Dolmetsch, K. Xu and R. S. Lewis, Nature 392 (1998) 933. 
[23] P. De Koninck and H. Schulman, Science 279 (1998) 227. 
[24] E. Oancea and T. Meyer, Cell 95 (1998) 307. 

[25] J. A. M. Borghans, G. Dupont and A. Goldbeter, Biophys. Chem. 66 
(1997) 25. 

[26] T. Haberichter, M. Marhl and R. Heinrich, Biophys. Chem. 90 (2001) 
17. 



13 



